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Abstract 

We develop here a simple yet versatile model for nuclear fragmentation in 
heavy ion collisions. The model allows us to calculate thermodynamic prop- 
erties such as phase transitions as well as the distribution of fragments at 
disassembly. In spite of its simplicity the model gives very good fit to recent 
data taken at the Michigan National Superconducting Cyclotron Laboratory. 
The model is an extension of a lattice gas model which itself has strong over- 
laps with percolation models which have been used in the past to compare 
with nuclear fragmentation data. 
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I. INTRODUCTION 



Nuclear fragmentation has been intensively studied for the purpose of measuring a pos- 
sible liquid-gas phase transition in nuclear matter [0,0. The most accessible quantity in 
experiments is a measurement of the yield Y{A) against A where A is the number of nucle- 
ons in the composite that emerges from the collision . A maximum in the fluctuation in 
the fragment sizes could be an evidence for critical phenomenon. In this paper we develop a 
schematic model for nuclear fragmentation. This is built on our previous work with a lattice 
gas model ^ which in turn was closely linked with much used percolation model for nuclear 
fragmentation [^,0. The model we propose can be related to phase transitions directly but 
we can also calculate the yield Y{A) against A. This yield against A can be calculated as 
a function of incident beam energy. There is just one parameter in the model, namely, the 
freeze-out density. This, however, is approximately known from past experiences in heavy 
ion collision theories. The model presented here involves only moderate computing and can 
be enlarged to include other effects. 



II. THERMODYNAMICS OF THE MODEL 

We consider a nearly central collision of approximately equal mass nuclei; because of 
nucleon-nucleon collisions the system reaches thermal equilibrium. We assume classical 
statistical mechanics is appropriate. In that case the canonical partition function of a n- 
particle system can be written in separable form: Z(can) = Zp(can)Zr(can) where Zp(can) 
is given by 

/n 
exp[-/5^p2/2m]d=^pi d^p.,. (2.1) 
1 

with f3 being the inverse temperature and pj being the momentum of particle i. The other 
part of the partition function is 

Zr{can) (X J exp[-/3^ v{rij)]d^ri d^Vn- (2.2) 

i<j 
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Here the f (ry ) is the potential between particles i and j. We approximate the configuration 
space part of the partition function Z^.(can) by the partition function of the lattice gas model 
I^J^. The motivation for introducing the lattice gas model is that it allows calculation of 
cluster distribution as in percolation model but also allows calculation of thermodynamic 
properties including liquid-gas phase transition in the same framework. 

We consider a simple cube lattice. Each lattice site may have one particle at the most. 
For a system of n particles the number of lattice sites, N, should be greater than or equal to 
n. When N = nwe have a nucleus in its normal volume. Thus the model considered here is 
limited to normal volume and larger. Cluster formation takes place at a freeze-out density 
which is expected to be about half of normal nuclear density so this is not a debilitating 
limitation of the model. 

We assume that the interaction is nearest neighbour type. This refiects the short range 
nature of nucleon-nucleon interaction. The number of nearest neighbours that a site has 
depends upon the dimension and the structure of a lattice and will be denoted by 7. For an 
infinite simple cube lattice we have 7 = 6. For finite systems that we will consider here the 
7 is less than 6 because of boundary effects. The energy of nearest neighbour interaction 
will be denoted by — e where e is positive and is related to the binding energy. 

We now obtain the equation of state for our system by first constructing the partition 
function. Let Nnn be the number of nn bonds in a specific lattice configuration; the energy 
carried by these bonds is then —eNnn- For a given N and n the partition function can be 
written as 

Z,,(can) = g{N, n, iV„„)e^^^"" (2.3) 
where g{N,n, Nnn) is a degeneracy factor satisfying 

Zr(can) was evaluated in in the Bragg- Williams approximation as well as in Bethe-Peierls 
approximation which is more accurate. It will be advantageous for us to rederive the Bragg- 
Williams result. In this approximation the number of nn bonds A^„„ is given and fixed 
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when and n are given. If one site is definitely occupied, the number of its 7 neighbours 
that are occupied on the average is 'yn/N. Since there are n nucleons in the system, the 
number of nn bonds is 'jn'^/2N where we have ensured that each bond is counted only once 
and assumed that both n and are so large that boundary effects can be neglected. The 
partition function is then 

Z,(can) = p^^^^^ (2.5) 

^ ^ {N-n)\n\ ^ ' 

The equation of state can be calculated by using P = kT [d In Z (can)/ dV)T = 

kT {d In Zr{caii)/dV) since Zp(can) has no V dependence. The volume V is given by 

V = a^N with = 1/po = 6.25 fm'^ where po is the normal nuclear density. The ground 

state volume of the nucleus consisting of n nucleons is Vq = a^n. Using Stirling's formula 

for A^!, n! and (A^ — n)\ it is easy to show that 

^ kT ^ N 1 ,n,. , , 

Using n/N = Vq/V = p/po we finally get 

^ = *rp„l„^-i<p„,(-^)^ (2.7) 

This equation of state has the same qualitative behaviour as the Van der Waals gas. For 
one mole of gas, the Van der Waals equation of state is 



NAkT a , , 



The lattice gas pressure goes to infinity as V approaches Vq. The Van der Waals gas 
pressure goes to infinity as V is squeezed to the value h. For large V both the equations of 
state approach the perfect gas limit. The critical point can be obtained analytically from 
equation (|2.7|). By setting dP/dp = d'^P/dp^ = at the critical point we obtain pc = 0.5po 
and kTc = je/A. Better calculations leave pc unchanged but reduces the value of Tc to 
1.1275e [0. An improved calculation for the partition function was performed in ^ where 
the equation of state of the lattice gas was compared with that for nuclear matter in a mean 
field Skyrme interaction. The qualitative behaviours are the same. 



III. CLUSTER DISTRIBUTION 



As in percolation model calculations we simulate an event by Monte-Carlo sampling and 
determine in each event the cluster distribution. There are two samplings to be done. We 
have to determine which of the sites the nucleons occupy and what their momenta are. 
The two samplings can be done independently of each other. In filling up the N sites with 
n particles the Boltzmann factor (eq. p.2[ ) needs to be taken into account. We do this 
in the following way. Starting with a empty lattice we put the first particle at random. 
Once this has been put in, the 7 sites which are immediate neighbours of this filled site are 
assigned a probability oc exp[/9e] whereas all other empty sites have probability proportional 
to unity. We now put the second particle according to this probability distribution. If at an 
intermediate step there are m empty sites we assign to each of these m sites a probability 
proportional to exp[g/5e] where q is the number of its nearest neighbours that are already 
occupied. The next filling is then done according to this probability distributiion. After the 
nucleons are all assigned their places among the N sites, the momentum of each particle is 
assigned according to Maxwell-Boltzmann distribution. 

Now let us look at the distribution of clusters. In a static lattice gas model in which the 
particles have no momenta, because of binding energy two neighbouring nucleons are always 
bonded and belong to the same cluster. The cluster is solely determined by this spatial 
configuration, as in the case of a percolation model. 

In our model, however, two neighbouring particles may not form a cluster because their 
relative kinetic energy might be higher than the attracting potential and the bond might be 
broken. It is natural to insist that two neighbouring nucleons are bonded if the following 
condition is satisfied. 

p^/2/i-e<0 (3.1) 

where p,- is the relative momentum and fi is the reduced mass. 

The frequency with which two nucleons appear in neighbouring sites is mostly deter- 
mined by the density. The probability of Pr/2/i exceeding the value e increases with tem- 
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perature since the momenta of individual nucleons are obtained from Monte-Carlo sampling 
of Maxwell-Boltzmann distribution at a given temperature. Hence, the probability that the 
two neighbouring nucleons are bonded decreases with increasing temperature, and the sys- 
tem becomes less compact at higher temperature as it should. A different parametrisation 
used in leads to similar effect. It is also clear that one important feature of pure per- 
colation models remains. Provided the density is not too small one will reach a percolation 
point in our schematic model at a certain temperature. This is the point when for an infinite 
system one infinite cluster just appears. For finite systems this signal for percolative phase 
transition appears with a second moment reaching a maximum which we discuss in the next 
section. 



IV. THERMAL AND PERCOLATIVE PHASE TRANSITIONS 

The power-law for emerging composites Y{A) oc A^'^ was first established by the Pur- 
due group in their experimental data. The Michigan experiment has identified 0] the 
minimum in r as a function of the bombarding energy in near central collisions of two ap- 
proximately equal ions. At this minimum the fluctuations are maximum indicating that 
this may be a critical point. In the framework of the model we have introduced this will 
coincide with a percolation point; does it also coincide with the thermal critical point? The 
two types of critical points need not overlap. 

A percolation point can be found from the behaviour of the second moment of cluster 
size defined as |llO[| : 



= (4.1) 

n 

where the P{s) is the average number of clusters of size s. The infinite cluster is excluded 
in the summation for an infinite system. At the percolation point 5*2 diverges in an infinite 
system and is at maximum in a finite system. Campi |^ was the first to exploit this feature 
by studying S2 at different numbers of fragments. 
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In our model (with cluster algorithm given in the density determines the temperature 
at which the percolation point is reached. In Fig. 1 we show 5*2 as a function of temperature 
for n=85 (to simulate ^°Ar+^^Sc collisions in ^) on lattices A^=5'^, 6^ and 7^ to simulate 
densities p = O.GSpo, 0.39po and 0.25po- In the figure the largest cluster is excluded in every 



event ||TT|. If we assume that the density is close to the critical density pc = 0.5po, then 5*2 
is maximum at the critical temperature Tc- The minimum in r is also observed at T = Tc 
and p = Pc- That is, the percolation point coincides with the thermal critical point at the 
critical density. Furthermore, it is found that when the density increases, the temperature 
at which the percolation point is reached increases and has the value of 1.55Tc at p = po. 



This is in full agreement with the results obtained from the Coniglio-Klein algorithm |T2 
that was derived mathematically by mapping the partition function into a percolation model 
|r3| . At this moment, it is not clear if there is any a priori reason for this remarkable 
result to happen. 

V. COMPARISON WITH DATA 

We want to compare our model with the data presented in where the data are fitted 



to a power-law and the deduced exponent r is plotted against the beam energy. In our model 
we can obtain the value of r as a function of temperature. To compare with experiments we 
need to relate the temperature to excitation energy which then can be related to the beam 
energy. This can be done to different degrees of sophistication the simplest of which is to 
assume completely classical limit. In that case there is no kinetic energy at zero temperature 
and the ground state energy per nucleon is —eN^^^/n where N'^^ is the maximal number 
of bonds possible. For a system of 85 nucleons the binding energy per nucleon is about 8.5 
MeV. Since is determined by the geometry it can be used to fix the value of e. At the 

temperature T the average energy per particle is 1.5kT — tNnn/n where A^rm, the average 
value of Nnni is obtained from the computer simulation. We can then write 

lkT + e{NZ^-Nnn)/n = e^ (5.1) 
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where e* is the excitation energy per nucleon. For equal mass non-relativistic nuclear colli- 
sions we have 



e* = ^beam/4 (5.2) 

where -Ebeam is the beam energy per nucleon in the laboratory (there is an implicit assumption 
here that all available energy is thermalized). Thus the temperature is related to the beam 
energy. Using this mapping we compare the calculated r (dashed line) with the experimental 
data in Fig. 2. In the calculation we used n = 85 and = 6^ for freeze-out density 
p = 0.39po- The fit is qualitatively correct but the r increases too fast with beam energy. 
This is clearly because we assumed full thermalization. A much better quantitative fit is 



obtained if we use the temperature extracted from experiments directly flgjlq] . In this 
approach, the tail of the proton spectrum in the laboratory is fitted by assuming that the 
proton has a Maxwell-Boltzmann distribution in a source which is moving in the laboratory. 



We take this mapping from the figure given in [|T^. The difference in mapping between 



temperature and beam energy in our pure theoretical approach (c.f. eqs. |5.1| and ^.2|) and 
the phenomenological approach is shown in Fig. 3. 

When this phenomenological mapping of temperature with beam energy is used the fit 
between the experimental data taken from P, P^JT7[ ] and our calculation is rather remark- 
able with one exception; the minimum value of r seen in experiment is lower than what we 
calculate. It has been pointed out that such low values of r probably signify that fragmen- 
tation took place in less compact geometry |18|. For example if clusters materialize from 



a toroidal or a bubble-shaped configuration the calculated value of r would be lower. The 
determination of such shapes however require a dynamical calculation. Quantum effect may 
also play an important role at low temperatures. 

Fig. 4 repeats the calulation of Fig. 2 except that we assumed at dissociation N = 7^ 
leading to p/po = 0.25. The fit is not as good as that for p/po = 0.39, especially the 
minimum in r is hardly discernible. Our results indicate that p/ po = 0.39 is a better choice. 
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VI. CONCLUSIONS 



The present model includes interactions between nucleons and the momenta of nucle- 
ons. It is more realistic than the percolation model which has been used in the past to 
look for signatures of phase transition. A significant feature is that the same model en- 
compasses thermal critical point and percolation points. The very good agreement between 
our calculations and the data taken at Michigan National Superconducting Cyclotron is 
very encouraging for our model. We feel that existing data are strongly suggestive of phase 
transition although more work needs to be done to establish this more firmly. This includes 
the measurements of various critical exponents and comparing them with calculations using 
models such as ours. Work in this direction has begun and is encouraging. 
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FIGURE CAPTIONS 



Fig. 1 The second moment 5*2 is plotted as a function of T /Tc at different densities. Here 
Tc = 1.1275e is the critical temperature in a lattice gas model. 

Fig. 2 The theoretical exponent r is compared with experimental data at different beam 
energies. The dashed curve is obtained by using the temperature calculated from eqs. 
( [5. ID and (|5.2|), and the solid curve is obtained by using the temperature fitted from 
experimental data. The solid circles are the corrected data taken from |]T^, the open 



circles are the uncorrected data taken from HI, and the crosses are taken from 17 



Fig. 3 The temperature is plotted as a function of beam energy. The dashed curve is 
calculated from eq. ( ^.1|) and (|5.2| ), and the solid curve is taken from where the 
temperature is obtained by fitting experimental data. 

Fig. 4 The same as the Fig. 2, but for p/po = 0.25. 
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